# This files estimates the discontinuity for variables expected to be unaffected to children's eligibilty
load("agg_day_year_all.rdata")
load("n_day_year_all.rdata")

## commented out for replication files 

#  install.packages("ggplot2")
#  install.packages("rmeta")
#  install.packages("reshape2")
#  install.packages("dplyr")
#  install.packages("rdrobust")

library(dplyr)
library(rmeta)
library(rdrobust)
library(ggplot2)
library(reshape2)

# write function for two samples proportion test

pr_mean_se <- function(var, by){
  p_1 <- var[by == 1]
  p_0 <- var[by == 0]
  mean_1 <- summary(p_1)[4] 
  mean_0 <- summary(p_0)[4]
  pr_mean <- mean_1 - mean_0
  pr_se   <- sqrt( mean_1*(1-mean_1) / length(p_1) + 
                     mean_0*(1-mean_0) / length(p_0) )
  p       <- 2*(1-pnorm(abs(pr_mean)/pr_se)) 
  n_under <- length(p_0)
  n_over  <- length(p_1)
  return(c(pr_mean, p, pr_se, mean_0, n_under, n_over))
}

# choose bandwidth based on when p < 0.05 parents with children 14 through 17 years or 19 through 22 years



plot_data <- 
  data_frame(est = rep(NA, 175),
             up  = rep(NA, 175),
             low = rep(NA, 175),
             year  = rep(c(2009, 2013, 2014, 2015, "Pooled"), 35),
             place = rep(1:5, 35),
             alpha = rep(c(rep(0.5, 4), 1), 35),
             var   = rep(rep(c("Parent \nvoted",
                               "Parent is \nnative",
                               "Parent is \nfemale",
                               "Parent has \nbachelor's",
                               "Child is \nfemale",
                               "Child is \noldest",
                               "Parent and child \ncohabit"), each = 5), 5),
             time  = rep(c("One day", "Two days" , "Three days", "One week", "One month"),
                         each = 35)) 

data_master <- 
  left_join(data_day_year, data_n_year, by = c("days", "year")) 


mock_vars <- c("par_vote",
               "par_native",
               "par_female",
               "par_degree",
               "chi_female",
               "chi_oldest",
               "par_cohabit")

day  <- c(1, 2, 3, 7, 30)
year <- c(2009, 2013, 2014, 2015)
for (j in 1:7){
  var <- mock_vars[j]  
  data <- 
    data_master %>%
    ungroup() %>%
    mutate(par_vote  = unlist(data_master[,var])) %>% 
    select(year, days, par_vote, n)
  for(i in 1:5){
    window <- day[i]
    data_close <- 
      data %>%
      filter(abs(days) < window) %>%
      mutate(voted     = round(par_vote*n),
             abstained = round((1-par_vote)*n),
             treated   = days > 0) %>% 
      select(days, treated, voted, abstained, year)
    
    data_voted <- 
      as.data.frame(lapply(data_close, 
                           function(x,p) rep(x,p), 
                           data_close[["voted"]])) %>%
      mutate(voted = 1)
    
    data_abstained <- 
      as.data.frame(lapply(data_close, 
                           function(x,p) rep(x,p), 
                           data_close[["abstained"]])) %>%
      mutate(voted = 0)
    
    data_close <- 
      rbind(data_voted, data_abstained)
    
    for(k in 1:4){
      est_year <- year[k]
      data_sub <- 
        data_close %>% 
        filter(year == est_year) 
      
      test <- with(data_sub, pr_mean_se(voted, treated))
      
      plot_data$est[(i - 1) * 35 + (j - 1) * 5 + k] <- test[1]
      plot_data$low[(i - 1) * 35 + (j - 1) * 5 + k] <- test[1] + qnorm(0.025) * test[3]   
      plot_data$up[ (i - 1) * 35 + (j - 1) * 5 + k] <- test[1] + qnorm(0.975) * test[3]   
    }
    
    test <- with(data_close, pr_mean_se(voted, treated))
    
    plot_data$est[(i - 1) * 35 + j * 5] <- test[1]
    plot_data$low[(i - 1) * 35 + j * 5] <- test[1] + qnorm(0.025) * test[3]   
    plot_data$up[ (i - 1) * 35 + j * 5] <- test[1] + qnorm(0.975) * test[3]       
    
  }
  print(j)
  
}

plot_data$time <- factor(plot_data$time, 
                         levels = c("One day", "Two days" , "Three days", "One week", "One month"))
plot_data$var <- factor(plot_data$var, 
                        levels = c("Parent \nvoted",
                                   "Parent is \nnative",
                                   "Parent is \nfemale",
                                   "Parent has \nbachelor's",
                                   "Child is \nfemale",
                                   "Child is \noldest",
                                   "Parent and child \ncohabit"))
plot <- 
  ggplot(data = plot_data,
         aes(x = 100*est, 
             xmin = 100 * low,
             xmax = 100 * up,
             y    = place)) + 
  geom_point(alpha = plot_data$alpha) + 
  geom_errorbarh(alpha = plot_data$alpha, height = 0) +
  facet_grid(var  ~ time, scales = "free_x") + 
  theme_classic() + 
  scale_y_continuous("", labels = c(year, "Pooled")) + 
  scale_x_continuous("Effect in percentage points") + 
  theme(panel.spacing = unit(0.5, "lines")) + 
  geom_vline(xintercept = 0, alpha = 0.3, linetype = 2)

